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Generating a unitary transformation in the shortest possible time is of practical importance to quantum in- 
formation processing because it helps to reduce decoherence effects and improve robustness to additive control 
field noise. Many analytical and numerical studies have identified the minimum time necessary to implement a 
variety of quantum gates on coupled-spin qubit systems. This work focuses on exploring the Pareto front that 
quantifies the trade-off between the competitive objectives of maximizing the gate fidelity T and minimizing the 
control time T. In order to identify the critical time T* , below which the target transformation is not reachable, 
as well as to determine the associated Pareto front, we introduce a numerical method of Pareto front tracking 
(PFT). We consider closed two- and multi-qubit systems with constant inter-qubit coupling strengths and each 
individual qubit controlled by a separate time-dependent external field. Our analysis demonstrates that unit 
fidelity (to a desired numerical accuracy) can be achieved at any T > T* in most cases. However, the optimiza- 
tion search effort rises superexponentially as T decreases and approaches T* . Furthermore, a small decrease 
in control time incurs a significant penalty in fidelity for T < T* , indicating that it is generally undesirable to 
' (— i ' operate below the critical time. We investigate the dependence of the critical time T* on the coupling strength 

Q | between qubits and the target gate transformation. Practical consequences of these findings for laboratory im- 

plementation of quantum gates are discussed. 
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I. INTRODUCTION 

The goal of controlling the dynamics of a quantum system in order to generate a target unitary transformation is 
both of fundamental interest and directly applicable to implementation of logic operations in quantum information 
processing [1]. Two strategies are commonly employed to design control fields that enact the desired evolution: 
(i) geometric techniques for analytically constructing control pulse sequences [2-A] and (ii) numerical methods 
employing optimal control theory (OCT) [5-8]. The effectiveness of the OCT approach has been demonstrated 
both for the idealized case of closed quantum systems undergoing unitary evolution and for open quantum systems 
whose dynamics are affected by coupling to the environment [5, 6]. Control fields producing quantum unitary 
transformations with a high fidelity have been successfully identified using OCT for a variety of models involving 
coupled-spin systems [9-17], molecular systems [18-20], and other physical realizations [21-28]. 

An important physical parameter for quantum computation is the control time T required to generate a target 
quantum gate. In general, decreasing the control duration helps to reduce the effect of decoherence resulting from 
the interaction of a quantum system with the environment [1]. Also, the gate error due to additive white noise in 
the control fields grows linearly with T [29], which means that shorter control times will enhance the robustness to 
this type of noise. Due to these considerations, the problem of identifying control fields that enact a target quantum 
gate to a specified fidelity in the minimum time, called time-optimal control (TOC) [2], is important for practical 
quantum information processing. The minimum time required to implement a multi-qubit gate is also related to 
the gate's complexity expressed as the number of one- and two-qubit gates necessary to construct the target unitary 
transformation [30]. TOC was originally formulated as a geometric problem of identifying the geodesic between 
two elements of the unitary group U(7V) and solved using Lie group methods and Pontryagin's minimum principle 
[2]. This analytical technique has been applied to identify control pulse sequences and associated values of the 
minimum time T* for generating quantum gates in two- and three-qubit NMR and other coupled-spin systems [2- 
4, 31, 32]. An alternative approach is to solve variational equations for the optimal time-dependent Hamiltonian 
H(t) under the constraint of finite energy using the quantum brachistochrone method [33, 34], which can provide 
improved T* values for some systems [34]. Other OCT-based algorithms utilize additional cost terms to penalize 
the control duration [35, 36]. In addition, several numerical studies have identified T* values approaching or even 
improving upon analytical results when employing OCT to design optimal fields [9, 10, 37]. 

In this work, we incorporate the goal of TOC by considering a more general problem of quantum Pareto op- 
timization [38] for the objectives of maximizing the gate fidelity T and minimizing the control time T. These 
two objectives compete with each other when T < T*, with the relationship between the best simultaneously 
achievable values of T and T constituting the Pareto front. Previously, Pareto optimization has been explored 
both theoretically and experimentally for the goal of discriminating between similar quantum systems [39, 40], as 
well as theoretically for maintaining persistent field-free control [41]. For TOC, some numerical simulations have 
sought the best fidelity value attainable at a given T [10, 37], but without explicitly investigating the Pareto front, 
especially in regions of high fidelity that are important for quantum information processing. We explore these 
Pareto fronts with the goal of identifying the relationship between the simultaneously achievable T and T values, 
as well as their dependence on the target unitary transformation and inter-qubit coupling strength. In order to 
numerically implement this analysis, we introduce the Pareto front tracking (PFT) algorithm, which sequentially 
(a) makes a small variation in one objective (here, decreases the control time T) and (b) searches for a control 
field that optimizes the second objective (here, maximizes the gate fidelity _F). Unlike the procedures described 
in Refs. [38, 39], the PFT method does not simultaneously optimize both control objectives, and thus may be less 
computationally expensive, especially when the value of one objective (here, T) is easily varied, but not easily 
optimized using OCT. The PFT algorithm is applicable to any such pair of objectives, for example, control field 
fluence and fidelity. In this work, we consider only the objectives of minimizing T and maximizing T . 

In addition to exploring the region of the Pareto front corresponding to T < T* (where the maximum attain- 
able fidelity is limited to values below 1), it is also of interest to understand how the optimization search effort 
(quantified as the number of algorithmic iterations needed to reach the optimum to a desired numerical accuracy) 
is affected when approaching the critical value T* from T > T* . For the objective of generating unitary trans- 
formations, the search effort was found to exhibit large variations with respect to the Hamiltonian structure [42]. 
In this work, we observe that the search effort to find optimal control solutions grows very rapidly as T decreases 
towards T*, even though a fidelity value arbitrarily close to the maximum T = 1 is, in principle, attainable at any 
T > T*. To facilitate the understanding of this behavior, we consider properties of the optimal control landscape, 
which is defined by the functional dependence of the physical objective (here, the gate fidelity J 7 ) on the applied 
controls [5, 6, 43—46]. For controllable quantum systems [47] with unconstrained control resources, the set of 
regular critical points on the unitary-transformation control landscape contains no local optima [46, 48, 49]. This 
property of the landscape topology is directly relevant to the optimization behavior, since local optima may act 
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as "traps" for a gradient-based search. The lack of traps on the control landscape has been verified with carefully 
conducted numerical simulations that ensured that no significant constraints were placed on the control fields [42]. 
Since the goal of TOC inherently involves limiting an important control resource (specifically, the control time T), 
it is possible that the favorable landscape topology may break down as the critical time T* is approached (fortu- 
nately, as we will show, this possibility does not materialize). Furthermore, it is of interest to explore how the local 
structure of the control landscape changes near T* and how this change is related to the rise of the search effort 
in this region. The PFT method introduced in this work is well-suited for exploring the landscape regions in the 
vicinity of the maximum fidelity while approaching T* from T > T* . In order to quantify how the local landscape 
structure changes upon approaching T*, we employ metrics similar to those developed in Refs. [42, 50, 51] and 
demonstrate their correlation with the search effort. 

The remainder of this paper is organized as follows. Section II presents the background and motivation for the 
current study, including the formulation of the optimal control problem, the optimization algorithm, metrics on the 
control landscape, model physical systems, the relationship between robustness to additive white control noise and 
control time, and the method for tracking the Pareto front. In Sec. Ill, we explore the fidelity-time Pareto fronts 
and the effect of control-time reduction on the search effort for two-qubit gates. The study is extended to three- 
and four-qubit gates in Sec. IV. Finally, Sec. V presents concluding remarks. 



H. BACKGROUND AND MOTIVATION 

A. Formulation of the control objective 

We consider an V-level closed quantum system whose evolution is governed by the time-dependent Schrodinger 
equation (in units where h = 1): 

i^^- = H({e k (t)})U(t,0), [7(0,0) = 1, (1) 

where H({sk(t)}) is the Hamiltonian, {ek(t)} are time-dependent external control fields, U(t,0) is the unitary 
propagator (time-evolution operator) from time t = to i, and 1 is the identity operator. We will use the shorthand 
notation U(t) = U (t, 0) for simplicity, where applicable. The propagator at some final time T is denoted as 
Ut = U(T) and is a functional of the control fields: Ut = UT({sk(t)}). We assume linear (dipole-type) coupling 
to the control fields: 

H({s k (t)}) = H + J2e k (t)H^ k \ (2) 

k 

where Hq is the field-free Hamiltonian and {H^} are the control-Hamiltonian operators. The quantum system is 
assumed to be evolution-operator controllable, which means that any desired Ut G U(iV) [or, Ut € SU(iV) for 
a traceless Hamiltonian] can be generated through the Schrodinger evolution (1) by some choice of control fields 
{sk{t)} at a sufficiently large time T [6, 8, 47]. The necessary and sufficient condition for evolution-operator 
controllability is that the operators {Hq, H^} generate the Lie algebra u(JV) [su(JV) for a traceless Hamiltonian] 
[47]. 

In circuit-model quantum computing, the goal is to generate specified unitary transformations that implement 
desired logic operations on a system of qubits. The corresponding control objective is to guide the system's final- 
time propagator Ut to match a specified unitary transformation W. A convenient mathematical formulation of this 
objective is to minimize the distance between Ut and W: 

V(U T ) = \\W- Ut\\hs = 2N- 2RcTr(IV t C/ T ), (3) 

where ||-X"||hs = Tr(X'X) is the squared Hilbert-Schmidt norm. The desired minimum T> = is achieved when 
Ut = W, and the maximum T> = AN corresponds to Ut = —W. It is often convenient to use the normalized 
distance V: 

^ (f/T) = \h v = \~ ^ ReTr ( wt ^), (4) 

so that T> takes values in the interval [0, 1]. A commonly employed gate fidelity T is related to the distance as 

T(U T ) = 1-V = ± + ^ReTr(W*U T ). (5) 
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The minimum distance (X> = T> = 0) at f/r = W corresponds to the maximum fidelity T = 1. The functional de- 
pendence of the objective on the control fields, i.e., T> = T>({ek{t)}) or, equivalently, T = J r ({efe(t)}) determines 
the optimal control landscape. 

OCT is often formulated by applying the variational principle to an objective functional, such as T (or, equiv- 
alently, T>), along with Lagrange multipliers to ensure satisfaction of the Schrodinger equation (1) as well as to 
impose a constraint on the control field fiuence [6, 1 8, 42] . An alternative approach [49, 5 1 , 52] is to consider small 
responses in the propagator U (t) due to changes in the control fields Se k , subject to the Schrodinger equation: 

d 

i-Q^U{t, 0) = H{t)SU{t, 0) + 5H(t)U(t, 0), (6) 
with the initial condition 5U(0, 0) = 0. Equation (6) can be integrated [49, 52] to give 

5U(t, 0) = -i [ dt'U(t,t')6H(t')U(t',0). (7) 



For the Hamiltonian of Eq. (2), the variation is given by 5H(t) = 2~^j. He 5Ek(t); by using this result in Eq. (7), 
one obtains [49] the variation of the propagator Ut with respect to the control £fc(t): 



SU 



T 



Se k (t) 



-iU T H^(t), (8) 



where H^ k \t) = W(t)H^ k) U{t). Combining Eqs. (3) and (8), we obtain the desired functional derivative of the 
distance V{Ut) with respect to £fc(i): 



SV 

-2RcTr 



5e k (t) 

= -2ImTr 



5e k {t) 
W ] U T H^ k \t) 



(9) 



The critical points of the control landscape (also referred to as extremal solutions) are control fields {sk(t)} that 
satisfy 

AT) 

0, Vfc and Vf G [0,T]. (10) 



<fefc(t) 

Quantum control landscape theory has shown [6, 46] that when (a) the system is controllable, (b) no significant 
constraints are placed on the control fields, and (c) the Jacobian in Eq. (8) is full-rank, dynamic critical points 
satisfying Eq. (10) occur only at kinematic critical points that satisfy VD(Ut) = 0. The values of the distance 
T> at the critical points are T> = 0, 4, 8, ... , AN [48, 49]. The optimality of a critical point can be determined by 
inspecting eigenvalues of the Hessian matrix H(t, t') = 5 2 V/5e(t)Se(t') [53]. Thus, the values T> = and T> = 
AN correspond to the global minimum (the Hessian is positive semidefinite) and global maximum (the Hessian 
is negative semidefinite), respectively. Moreover, it has been shown through analysis [48, 49] and numerical 
simulations [42] that, under the conditions (a)-(c) above, all the intermediate critical points (i.e., T> = 4, . . . , AN — 
A) have a saddle-point topology (the Hessian has positive, negative, and zero eigenvalues), meaning that no local 
maxima or minima exist on the landscape. However, when control resources are severely constrained (e.g., by 
limiting the control time T as considered here), the trap-free landscape topology is no longer guaranteed. 

Values of the distances and fidelity of Eqs. (3), (4), and (5) depend on the global phase of the transformation 
Ut- Since this global phase is physically irrelevant for a given quantum gate, a phase-independent version of the 
distance (or fidelity) is often employed instead [49]. In particular, a normalized phase-independent distance can be 
defined as [54] 

S( u t) = 4?min2?(e^l7 r ). (11) 
The minimization over the global phase in Eq. (11) can be easily carried out to obtain: 

g(U T ) = l-^\Tr(W^U T )\. (12) 
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Correspondingly, the minimum value Q = is attained when Ut = e"^W for any phase <j>. Note that for a traceless 
Hamiltonian, when both W and Ut must be in SU(iV), the phase <fi can take only discrete values corresponding 
to solutions of the equation e lN( ^ = 1. The topology of the control landscape for the distance Q of Eq. (12) is 
very similar to that for the distance T> of Eq. (3), and an analytical formula for the functional derivative 8Q / 5e k (i) 
can be obtained analogous to Eq. (9) [49]. Optimization of the phase-independent distance Q will be considered 
in Sec. IIIB, where we study how searches with different initial control fields converge to optimal solutions corre- 
sponding to different values of the global phase. Optimization of the phase-dependent distance T> is considered in 
the remainder of this paper (of course, minimizing the distance T> is equivalent to maximizing the fidelity F). 



B. Optimization procedure for control of unitary transformations 

A variety of deterministic first-order algorithms, including the Krotov method [ 1 8, 37, 55] , GRAPE algorithm [9, 
10], and D-MORPH (diffeomorphic modulation under observable -response-preserving homotopy) [42, 56] have 
been employed for optimization in control of quantum unitary transformations. A recent work [42] demonstrated 
that these algorithms share a common fixed point topology and common bounds on their convergence rates. Also, 
methods for comparison and benchmarking of various quantum control algorithms were presented in Ref. [57]. 
Second-order algorithms such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton method [57, 58] 
and the Newton-Raphson method [59] were also utilized recently in quantum optimal control. In this work, we 
employ the D-MORPH method. 

In D-MORPH, a variable s (referred to as the algorithmic index) is introduced to label the progression of the 
optimization, rendering the control fields: {sk(t)} —> {s k (s, t)}. The objective value T> depends on s through its 
functional dependence on the set of control fields, i.e., T>(s) = T>({s k (s, t)}). Thus, the change in the objective 
value corresponding to a differential change ds is given by cTD = (dT>/ds) ds, where 

^=f T dtY 5V d£kM (13) 

As the goal is to minimize T>, we require that dV/ds < 0, which is guaranteed when each control e k (s, t) satisfies 
the differential equation 

de k (s,t) _ 5V 
ds 8e k (s,t)' 

The functional derivative on the right-hand side of Eq. (14) can be evaluated using Eq. (9). In numerical simula- 
tions, we determine control fields at each iteration by solving Eq. (14) using a fourth-order Runge-Kutta integrator 
with a variable step size incorporated into MATLAB (routine ode4 5) [60]. The initial set of fields {e k (0,t)} is 
selected randomly for each optimization run, as described in Sec. II D below. In the simulations, we evaluate the 
normalized distance T> and specify the convergence threshold for the optimization. Specifically, the D-MORPH 
procedure described above is performed until either (a) the desired convergence criterion T> < 10~ 8 is reached or 
(b) the improvement of the distance value satisfies \D{s + ds) — 2?(s)| < 10~ 6 X>(s), indicating that an extremal 
value of T> has been reached. Sufficient algorithmic iterations are allowed to prevent premature termination before 
reaching an extremal value of T>. 



C. Metrics of landscape structure 



The issue of the search effort growth as T — > T* is important because the computational cost of identifying 
optimal control fields near or at the critical time may become prohibitively high, particularly as the number of 
controlled qubits increases. We are also interested in exploring the relationship between the number of D-MORPH 
algorithmic iterations needed to converge to an optimal solution and metrics that quantify local properties of the 
control landscape. One such metric is the path length of the search trajectory for a D-MORPH optimization. The 
search starts out from an initial set of control fields {efc(0, t)} at the algorithmic index value s = and progresses 
in steps s — s- s + ds until the trajectory ends at a set of optimal control fields {e£} = {e k (s*, t)} with s = s*. The 
path length A(s) of the search trajectory from s = to s is defined as 



dt 



~ de k (s',ty 
ds' 



1/2 



(15) 
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The total path length traveled along the search trajectory to reach the optimum is A* = A(s*). In the absence of 
local traps, a D-MORPH search monotonically converges towards the optimum, resulting in a one-to-one mapping 
between the algorithmic index value s and the objective value T> for a given search trajectory. Therefore, it is 
possible to cast the path length as a function of the distance: A = A(T>). The value of A (2?) quantifies how much 
the control fields have to change from the initial guess {ek{0, t)} to achieve the objective value T>. 

In Ref. [42], it was found that metrics of the local landscape structure near the optimum can be useful for 
predicting the search effort. In particular, one can employ the slope metric £(s), which is given by the L2-norm of 
the landscape gradient V e 2? evaluated at an algorithmic index value s: 



E( S ) = ||V £ I?( S )|| 2 = |^^ 



dt 



5V 



5e k {s,t) 



1/2 



(16) 



During the D-MORPH optimization, the derivative of the fcth control field with respect to s and the landscape 
gradient along this field are related by Eq. (14). Substituting Eq. (14) into Eq. (15) and using Eq. (16), we obtain 
[61]: 



AGO = 4= [ S ds'X(s>) & E(a) = Vf 
vT Jo 



c ZA(g) 
ds ' 



(17) 



This simple relationship shows that the slope metric S(s) is proportional to the rate of change of the path length 
A(s) along the D-MORPH search trajectory. As was the case for optimization searches in Ref. [42], we will show 
that the decrease of the slope metric S(s) (i.e., the increase of the landscape "flatness") at an algorithmic iteration 
close to the landscape optimum (e.g., for V w 1CP 6 ) con-elates with the increase of the search effort. 



D. Model control systems 



Various types of coupled-spin systems have been considered in TOC studies, with many works using models 
relevant to the liquid-state NMR, where spins are coupled via an Ising-type interaction and coupling strengths 
are much smaller than differences between spin frequencies [3, 9, 10, 34]. Other models such as anisotropic 
controllable inter-qubit couplings [33, 37] have also been studied. In NMR models, control fields typically address 
each spin separately with independent x and y polarizations [9, 10]. 

In this work, we study TOC of generic coupled-spin model systems motivated by implementations of quantum 
computing in physical devices such as semiconductor quantum dots [62]. In our model, each qubit has a character- 
istic transition frequency and is controlled by a separate field with only the x polarization; exchange interactions 
between qubits are of the Heisenberg type with fixed isotropic coupling strengths. All system and control parame- 
ters are expressed in dimensionless units. 

Generally, we consider a system of n qubits (the Hilbert space dimension is N — 2"), with the model Hamilto- 
nian of the form (2). The field-free Hamiltonian Hq is given by 

n n—X n 

H = J2"kSi k) + J2 E J^SW-SW. (18) 

k=l k=l j=k+l 

Here, the operator Sa (a = x, y, z) denotes the tensor product of the spin operator for the fcth qubit with identity 
operators for all other qubits: 

S ( a k) = 1 2 ® • ■ ■ ® t 2 <g> S a I2 ® ■ ■ ■ ® I2, (19) 

V v ' * ' 

fc— 1 n—k 

where the spin operators are S = (S x , S y , S z ) = \(<j Xl a yi a z ), in terms of the Pauli matrices, and II2 is the 
2x2 identity matrix. Each qubit has a unique transition frequency 0;^ (corresponding to the presence of a static 
magnetic field in the z direction in the spin model), and isotropic coupling strengths J^ k j> between pairs of qubits 
are constant. In the simulations reported here, we used model systems with up to four qubits with frequencies 
LOk = 20, 24, 30, 40 and coupling constants j( fe > J ) in the range from 0.08 to 400. This broad range of coupling 
strengths represents the freedom inherent in considering coupled-spin systems in contexts other than NMR, such 
as semiconductor quantum dots, where interactions between qubits may be tuned by application of electric fields 
[15, 62-64]. 



7 



The control Hamiltonian H c (t) corresponds to the application of a separate time-dependent control field polar- 
ized in the x direction to each individual qubit: 

n n 

H c (t) = $>(^ =£e*(*)4 fc) , (20) 
fe=i fe=i 

where £fc(t) is the control field applied to the kth qubit and the operator is defined by Eq. (19). In the 
optimization procedure, each control field is labeled by the algorithmic index s (see Sec. II B above). The fluence 
of the fcth field at the index s is given by 

f k (s) = [ dtel(s,t). (21) 
Jo 

At the start of the optimization (s = 0), each field is initialized in the parameterized form: 

M 

e k (0,t)=A(t)J2sm(r] i t + (p i ), (22) 
i=i 

t G [0, T]. Here, A(t) = Aq exp[— 8"7r(t — T/2) 2 /T 2 ] is the Gaussian envelope function, frequencies {rji} (corre- 
sponding to M spectral components of the field; we usually use M = 10) are randomly selected from a uniform 
distribution on [0, O] (with SI being the largest transition frequency in Ho), {(fii} are random phases on [0, 2ir], The 
normalization constant Aq is chosen so that the fluence /fc(0) of each initial field is equal to 1. 

The parameterized form (22) is used only for initial control fields. At each step of the optimization algorithm af- 
ter the initialization (i.e., for s > 0), the value of the kth control field at each point on the time discretization mesh is 
allowed to vary freely and independently. This flexible set of control 'knobs' {Ek{s, At), £fc(s, 2 At), . . . , £fc(s, T)} 
allows the fluence of the kth field to vary freely when s > 0. The time step At is chosen such that At < 7r/(2fi), 
corresponding to the Nyquist frequency cjn = 7r/At > 2ft. In agreement with the Nyquist-Shannon sampling the- 
orem, this criterion was found to be sufficient to ensure that the time discretization does not affect the reachability 
of the global optimum [42, 51]. 



E. Relationship between control time and robustness to additive white control noise 

An important motivation for finding the minimum time necessary to enact a target unitary transformation is to 
improve the robustness of the gate operation to noise in control fields. Specifically, additive white noise (AWN) 
in optimal control fields induces gate errors that are linearly proportional to T [29]; a brief outline of this analysis 
is provided here. In the presence of additive noise, the actual control field is given by e(t) + £(t), where £(t) 
is a classical stochastic variable. For white noise, £(t) has zero mean: E{£(t)} = 0, and is delta-correlated: 
E{£(t)£(t')} = <7 2 <5(t — t'). Here, E{-} denotes the statistical expectation value over all noise realizations and a 2 
is the variance of the noise amplitude distribution. For AWN in multiple control fields, in general, we should also 
consider cross-correlations: E{£fc(f)£j(t')} = <r 2 f3kjS(t — t'), where < /3kj < 1 and fikk = 1 (if noise processes 
in different control fields are independent, then /3k j = 8kj)- 

According to the analysis in Ref. [29], for weak AWN in optimal control fields {e%(t}}, the statistical expectation 
value of the normalized distance T> is approximated (by expanding up to the second order in the noise amplitude) 
as 

ni>} * \^ E fa f dt % (*> *)• < 23 > 

k,j=i J ° 

Here, H*(t,t') denotes the Hessian matrix of T> evaluated at the optimum, and the diagonal elements (i.e., for 
t = t') of its blocks are time-independent [29, 49]: 

LJ* (f f] ^ 
Hfe ' (M) - Se k (t)S Sj (t) 



2N~ 



-Tr 



(24) 
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Substituting Eq. (24) into Eq. (23), one obtains an expression that reveals the linear dependence of the expected 
gate error on T: 



^ w jN a2T E ^ Tr 



4iV 



(25) 



fej=i 



For the control Hamiltonian of Eq. (20) employed here, Tr[H^ k) H^ j) } = Tr[Si k) Si j) ] = (N/4)5 kj . Using this 
result (together with /3kk = 1) in Eq. (25), we finally obtain: 



E{V} 



1 

16' 



a 2 nT. 



(26) 



Equation (26) shows that the error in the objective value is expected to grow linearly both in the number of qubits 
n and control time T. Thus, it is practically important to minimize T when AWN in control fields is present. 



F. PFT algorithm for quantifying the trade-off between minimizations of distance and control time 

The PFT algorithm introduced here is designed to identify the Pareto front for the dual objectives of minimizing 
T> and T, as well as to explore the corresponding domain of the optimal control landscape. We use PFT to move 
along the Pareto front by identifying optimal control solutions corresponding to different values of T. Analogous 
to the algorithmic index s describing the progression of a D-MORPH search, we define the indexing variable p to 
describe the progress of the PFT algorithm. The PFT algorithm works as follows: 

1. Select a starting value of the control time, Tq = T(p = 0). Then run the D-MORPH optimization (as 
described in Sec. II B), starting from a set of randomly selected initial control fields {Ek{p = 0, s = 0, t)}, 
until it converges to a set of optimal fields {sk(p = 0, s = s*, t)} that minimizes T>. 

2. Reduce the value of the control time T, so that T(p + 1) = T(p) — AT, where AT is an increment on the 
order of AT < 0.01T. 

3. Resample each of the optimal control fields in the set {ek{p, s*,t)} on the updated time interval [0, T(p + 
1)]. Employ the resulting set of fields as the initial guess {Ek(p + 1, s = 0, t)} for the next D-MORPH 
optimization that proceeds to identify the next set of optimal fields {sk(p + 1, s = s*,t)}. 

Steps 2 and 3 are repeated until the D-MORPH optimization can no longer attain the desired value of T>. Specif- 
ically, we set the PFT "stop value" to T> = 1CP 2 when the goal is to explore the "competitive" part of the Pareto 
front, and to T> — 1CP 8 when the goal is to only determine the critical time T*. The process of decreasing T in 
small increments coupled with the use of optimal control fields as the initial guess for the D-MORPH run in the 
next PFT iteration biases towards identifying families of related control solutions along the Pareto front. Therefore, 
running multiple PFT trajectories beginning from different random initial control fields at various values of To may 
be useful for proper identification of the Pareto front. 



m. EXPLORING THE PARETO FRONT FOR CONTROL OF TWO-QUBIT GATES 

In this section, we explore the Pareto front for fidelity- and time-optimal control of unitary transformations 
in two-qubit systems. In Sec. Ill A, we study in detail how the optimization search effort, landscape metrics, 
and optimal control fields change along the Pareto front for the controlled NOT (CNOT) target gate and one 
representative set of system parameters. In Sees. MB and IIIC, we identify Pareto fronts and determine critical 
times for a variety of target gates and inter-qubit coupling strengths. 



A. Optimization search effort and Pareto front exploration for CNOT gate 

As an illustrative case on which to conduct a detailed examination of the search effort dependence on the control 
time and properties of the distance-time Pareto front, we consider the objective of performing the CNOT gate in 
the two-qubit model system (n = 2) with u>i = 20, 0J2 = 24, and j' 1 ' 2 ' = 0.8. Since the system Hamiltonian is 
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traceless, the target transformation W is defined as the CNOT gate with a global phase factor chosen so that W is 
in SU(4): 



ir< 
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FIG. 1. (Color online) The search effort and metrics of the control landscape as functions of the control time T for D-MORPH 
optimizations performed at select fixed values of T. The target gate is Wcnot- (a) Search effort (black circles with solid line, 
left-side ordinate) and slope metric (the Z/2-norm of the landscape gradient) E(s) evaluated at the s value corresponding to 
V w 10 -6 (red triangles with dashed line, right-side ordinate), (b) Total fluence /* of optimal control fields (black triangles 
with solid line, left-side ordinate) and total path length A* along the search trajectory (red circles with dashed line, right-side 
ordinate). Circles and triangles show average values and error bars denote the left and right standard deviation over the sample 
of 10 D-MORPH searches with randomly selected initial control fields performed for each value of T. 



1. Search effort of D-MORPH optimizations and its relationship to landscape metrics 

For the selected target gate and control system, the critical time value T* « 4.12 was found by running PFT 
trajectories (see Sec. Ill A 2 below for details). In order to study how the control time affects the optimization 
search effort, D-MORPH trajectories were obtained for 15 values of T between T — 40 and T = 4.12. For each 
value of T, 10 D-MORPH optimization runs beginning from different random initial fields of unit fluence were 
performed. All optimizations reached the desired objective value T> < 10~ 8 and no trapping or slowdown of 
searches at a suboptimal distance value was observed (including runs for T = T*). The search effort, however, 
increased as T was made smaller, especially for T < 6. This behavior is shown in Fig. 1(a) as a plot of the 
number of D-MORPH algorithmic iterations (averaged over 10 searches started from random initial fields) versus 
T. In particular, for T < 6, the search effort increases superexponentially (note the logarithmic scale of the left- 
side ordinate) and is well approximated as cxp(aT~ b + c), where a w 2.1 x 10 7 , b ~ 11.0, and c w 6.1. The 
corresponding decrease in the slope metric E(s) evaluated at the s value corresponding to V w 10~ 6 is also shown 
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in Fig. 1(a) (with values on the right-side ordinate). The complementary trends for the search effort and gradient 
norm indicate that the control landscape in regions near the optimum becomes "flatter" as T decreases. 

To examine how the distance between the initial and optimal fields grows as T decreases towards T*, we 
consider the total fluence of optimal control fields, /* = J^k /fc( s *)> and the total path length A* = A(s*) traveled 
along the search trajectory to reach the optimum. These quantities (averaged over 10 searches started from random 
initial fields) are plotted versus T in Fig. 1(b). Both /* and A* rise as T decreases, and this rise significantly 
accelerates for T < 8. Since all D-MORPH searches began from unit-fiuence fields, the value of /* is an indicator 
of the distance between the initial and optimal fields, which explains the similarity in the trends of /* and A*. The 
behavior of the fluence as a function of the control time can be qualitatively explained by a simple example of a 
two-level system driven by a resonant control field. The rotation angle produced by the control Hamiltonian on 
the Bloch sphere during time T is CIrT, where flu is the Rabi frequency which is proportional to the control-field 
amplitude e - If the goal is to generate the same rotation as T changes, the optimal-field amplitude should scale as 
£g oc 1/T. Since the fluence can be approximated as / ~ \e\T, the optimal-field fluence scales as /* oc 1/T. Of 
course, for a multi-qubit system controlled by several external fields, the dynamics is much more complicated (in 
fact, even for a single qubit, 1/T scaling is periodically modulated by the effect of free evolution). Nevertheless, 
this simple picture helps to explain qualitatively why the fluence rises when T decreases. The non-monotonicity 
of /* and A* as functions of T, seen in Fig. 1(b), is explained by the fact that, as T changes, free evolution takes 
the system closer to or further away from the target, thus modulating the amount of control-field energy required 
to reach the optimum. 




FIG. 2. (Color online) Path length difference A(s*) — A(s) along the optimization trajectory as a function of the objective T>. As 
the optimization progresses, the path length difference decreases towards zero along with T>. The target gate is Wcnot- Each 
value of the control time T is denoted by line style and color in the legend. The order of trajectories in the legend corresponds 
to that on the figure, with the T* trajectories (green, solid line) at the top, and T increasing from top to bottom. Different 
trajectories for the same value of T correspond to 10 D-MORPH searches with randomly selected initial control fields. 



Further insight into the dependence of the control landscape structure on the control time can be gained by exam- 
ining the path length A(s) accumulated as the D-MORPH search progresses. The value of A increases monotoni- 
cally along the search trajectory from s = to s = s*; correspondingly, the path length difference A(s*) — A(s) 
indicates the extent to which the control field has reached its optimal form. Figure 2 shows A(s*) — A(s) as a 
function of the objective value T> for optimization trajectories with selected values of T, with 10 trajectories cor- 
responding to different random initial fields shown for each T. This plot illustrates the dependence of the search 
trajectory on T. For all searches with T > T*, the path length difference follows a similar power law (which 
appears as a linear change on the log-log plot) over the range from T> ~ 10~ 3 to T> ~ 10~ 7 . In contrast, the 
searches with T = T* approach the final path length A(s*) much more slowly until T> becomes very close to the 
optimum; there, the path length changes very quickly, as A(s*) — A(s) drops by three orders of magnitude between 
T> ~ 2 x 10~ 8 and T> ~ 10~ 8 . This result shows that at T*, large changes of the control field occur very near the 
optimum. 
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FIG. 3. (Color online) The normalized distance T> plotted versus the control time T for 10 separate PFT trajectories denoted 
by shape and color. The target gate is M^cnot- The PFT trajectories were started at different values of To (given in the legend) 
with different randomly selected initial control fields, but are very closely aligned, which suggests the existence of a unique 
Pareto front. 



2. PFT results 

The PFT procedure described in Sec. II F was performed for the purposes of (a) identifying the critical time T* 
(determined numerically as the minimum control time at which T> < 10~ 8 can be achieved) and (b) following the 
Pareto front to obtain the best attainable T> value as a function of T in the "competitive" region of T < T*. A 
total of 10 PFT trajectories beginning from random initial control fields at different values of To (ranging from 
Tq = 4.2 through T = 6) were generated for the target gate Wcnot and the same two-qubit system as used in 
the D-MORPH optimizations in Sec. Ill A 1 above. Four of these trajectories were run along the Pareto front by 
decreasing T until it was impossible to attain T> < 10~ 2 , while the remaining trajectories were run only to values 
of T ranging from T = 3.8 to T = 4 due to computational expense. All 10 trajectories, shown in Fig. 3, are very 
closely aligned for all values of T, with the critical time value estimated as T* = 4. 12 ±0.01 (i.e., the dispersion of 
the T* value is on the order of the PFT step AT). This near coincidence of different PFT trajectories indicates that, 
for a given quantum system and a target gate, there exists a unique distance-time Pareto front. Correspondingly, 
the minimum control time needed to enact a target unitary transformation with a high fidelity appears to be an 
inherent property of the controlled quantum system, independent of the path taken to identify T*. 

The Pareto front for T < T* has the unfavorable property of an extremely steep slope, i.e., a small decrease in 
T below T* results in a very large increase in T>. Specifically, T> rises more than three orders of magnitude from 
T> < 10~ 8 to T> > 10~ 5 with a relatively small decrease in T from T = 4.12 to T = 4.0. Beyond this steep region, 
the distance growth moderates with decreasing T, such that T> ps 10~ 2 can still be obtained at T = 3.0. However, 
since a fault-tolerant quantum computation requires very low gate error rates (typically, less than 10~ 4 ) [1, 65], the 
steep slope of the Pareto front immediately below T* presents a fundamental limitation on gate implementation 
times. Furthermore, uncertainty introduced under experimental conditions would make operation near T* difficult 
because small errors in control time could cause substantial decreases in attainable fidelity. 

In addition to being a reliable numerical method to identify the distance-time Pareto front (and, in particular, de- 
termine the value of T*), the PFT algorithm generates families of related control fields that minimize the objective 
T> for control times within a selected interval. While all PFT trajectories achieve essentially the same minimum 
value of 2? at a given T, each trajectory produces a distinct family of optimal control solutions. The "evolution" of 
optimal control fields within such a family, corresponding to the change of T along a PFT trajectory, is visualized 
in Fig. 4. Specifically, two families of optimal fields obtained for PFT trajectories initialized at T = 4.4 and 
T = 4.6 are shown in Figs. 4(a), (b) and 4(c), (d), respectively. For each of these two families, pairs of fields 
(corresponding to separate control fields acting on two qubits) are shown for several values of T ranging from 
T = T to T = 3.0. 

The optimal control fields obtained for different PFT trajectories have distinct shapes, which is expected based 
on the existence of an infinite number of optimal solutions [66] . Nevertheless, the fields share a number of common 
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FIG. 4. (Color online) Optimal control fields {el(t)} for the target gate VKcnot- The control fields for the first qubit (k = 1) 
are shown in (a) and (c), and the control fields for the second qubit (k = 2) in (b) and (d). These plots show fields obtained 
for various values of the control time T, ranging from T = To (the starting point of the PFT trajectory) to T = 3 (where the 
best attainable objective value is T> ss 0.01). Each value of T is denoted by color in the legend, with yellow (lightest gray) 
corresponding to T = T* = 4.12. The shown sets of fields represent two PFT trajectories started at To = 4.4 ((a) and (b)) 
and To = 4.6 ((c) and (d)); these PFT trajectories are presented among several others in Fig. 3. An apparent increase in the DC 
field component (i.e., a positive shift of the entire field) as T decreases is an artifact of the perspective of the three-dimensional 
plot; the field power spectra (not shown) do not reveal a significant zero-frequency component. 



features. First, field amplitudes increase as T decreases. As explained in Sec. Ill A 1, this amplitude behavior is 
needed to maintain the required rotation angle as T changes, and the associated change of the field fluence roughly 
follows 1/T scaling. Second, for T < T*, each field exhibits spikes at t = and t = T that grow in amplitude 
as T decreases. For T < 3.9, these spikes are the largest amplitude features of the control fields. Besides the 
two families shown in Fig. 4, this field feature was observed for all other PFT trajectories that were run below T*. 
Similar characteristics of optimal fields were also observed for TOC of a three-qubit system in Ref. [10]. We made 
no attempt to impose any control constraints that would suppress such spiky features; the extreme difficulty of 
the distance minimization in the region T < T* indicates that further constraints could adversely affect attainable 
values of T> and thus preclude reaching the genuine distance-time Pareto front. 

The results presented in this section suggest that the location of the Pareto front is essentially independent of the 
PFT trajectory taken for the present numerical model. Thus, for other selections of the target gate and/or system 
parameters considered in Sees. IIIB and IIIC below, Pareto fronts were identified by running only one complete 
(i.e., followed until T> < 10 -2 becomes unattainable) PFT trajectory. For cases where only identification of T* 
was desired, the PFT algorithm was stopped soon after T> < 10 -8 was no longer attainable, and at least three 
trajectories started with different initial random fields were run in this fashion in order to verify the obtained value 
ofT*. 
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FIG. 5. (Color online) The search effort versus the control time T for five PFT trajectories (denoted by distinct symbol shapes 
and colors) and D-MORPH searches with random initial fields (black squares with error bars denoting the average and standard 
deviation, respectively, over the set of 10 runs). At each value of T < 4.45, the PFT search effort is sizably smaller than the 
average D-MORPH effort. 



3. Efficiency of the PFT algorithm 

In Sec. IIIA1, a superexponential increase of the optimization search effort (in terms of algorithmic itera- 
tions) was observed as the control time T decreases and approaches T* . The D-MORPH searches considered in 
Sec. Ill A 1 were initialized at randomly selected fields for all values of T, with typical initial objective values 
T>(s = 0) ~ 0.5. In contrast, the PFT algorithm employs random initial fields only for the search with the starting 
value of the control time, To = T(p = 0); every consequent search with T(p + 1) < Tq is initialized at the fields 
{ek{p, s*, t)} that are optimal for the preceding search with T(p) = T(p+1) + AT. Correspondingly, D-MORPH 
searches along a PFT trajectory begin, forp > 0, with initial objective values T>(s = 0) ~ 0.01 to 0.05 and thus 
may be expected to reach V < 10~ 8 with a smaller number of algorithmic iterations. Here, we investigate to what 
degree the PFT algorithm can lower the search effort, as compared to D-MORPH optimizations with random initial 
fields. 

To directly compare the two methods, we ran (a) five PFT trajectories from To = 4.6 to T* = 4.12 and (b) 
sets of 10 D-MORPH searches with random initial fields at selected values of T E [4.12, 4.6]. Figure 5 presents 
the search effort as a function of T for each PFT trajectory (colored circles, diamonds, triangles, crosses, and 
x's) and for each set of D-MORPH searches (black squares and error bars indicating the average and standard 
deviation, respectively, over the set of 10 runs). The results show that, while different PFT trajectories exhibit 
varying convergence speeds, at each value of T < 4.45, the search effort for PFT is significantly lower than that 
for randomly initialized D-MORPH, typically by a factor of 2 to 4. At T* = 4.12, the convergence of the best PFT 
trajectory was faster than that of the average D-MORPH search by a factor of ~ 8; even the worst PFT trajectory 
outperformed the average D-MORPH search by a factor of ~ 2. Thus, the search effort can be substantially 
lowered if a D-MORPH search is initialized at control fields that are optimal for a related search (here, with a 
slightly larger value of T) rather than at random fields. 

It is also of interest to compare the PFT algorithm to methods that aim at minimizing T by including penalty 
terms into the cost functional [35, 36]. For the latter approach, a recent work [36] assigned a priori weights to 
a term in the cost functional that penalizes the control field fluence and control time. It is known [38] that OCT 
algorithms employing such cost functionals are typically incapable of identifying the genuine Pareto front. Indeed, 
the algorithm used in Ref. [36] converges to a value of T that overestimates the actual critical time, while also 
underestimating the achievable gate fidelity. A PFT trajectory that we ran for the quantum system and target gate 
used in Ref. [36] identified T* « 1.80 with J" > 1 - 10~ 8 , while Ref. [36] reported convergence to T « 2.01 with 
Jk1-7x 10~ 6 . These results suggest that the PFT algorithm is better suitable for accurately determining the 
fidelity-time Pareto front and, in particular, the true value of T* than methods that employ time-penalizing terms 
in the cost functional. 
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B. Pareto fronts and critical times for various target gates 



Simulations presented in this section examine how the distance-time Pareto front and, in particular, the critical 
time are affected by the choice of the target transformation W for the same two-qubit system as considered in 
Sec. Ill A. Pareto fronts and values of T* for various target gates are identified using the PFT procedure, as 
described in detail above. In addition to the Wqnot gate of Eq. (27), we considered a number of gates, all 
incorporating appropriate global phase factors to make them elements of SU(iV). The SWAP gate has the form 
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We also consider the square root of SWAP (an entangling two-qubit gate), given by 
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Matrix elements of the Quantum Fourier Transform (QFT) gate generalized for n qubits are given by 



W QF T,nU> k ) = -J^ e 
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e 2m/N ( reca n that N = 2 »). For two-qubit systems (n = 2), the QFT gate 
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We also consider a similar gate denoted as W-qft',?^ which is given by Eq. (30) with j, k = 1,2, 
we consider controlled phase (CPHASE) gates, of the form 



Wcphase(u) = e 
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. , N. Finally, 



(32) 



with the phases a = tt and a = tt/2. 

In addition to the explicit forms of the CNOT, SWAP, VSWAP, QFT, QFT', and CPHASE gates presented 
above, we also consider the effect of changing the global phase of the target gate. In earlier work by Schulte- 
Herbriiggen et al. [10], different values of T* were found for distinct global phases of the three-qubit QFT gate 
in a linear spin-chain system with Ising-type interactions. Identifying critical times that correspond to all possible 
values of the target gate's global phase is important in the context of TOC. For an n-qubit gate W E STJ(AT), 
transformations with N distinct global phase values are allowed: 



W{<f 3m )=e^W 1 

<j) m = 2mir/N, m = 0, 1, . 



.,N-1. 



(33) 



For two-qubit gates considered in this section, four distinct global phase values are allowed in SU(4): <\> m = ran/ 2, 
m = 0, 1, 2, 3. The critical time for the transformation W(<fi m ) is denoted as T*((f> m ). 

The effect of global phase on the location of the distance-time Pareto front is demonstrated in Fig. 6 for the 
SWAP and QFT gates with = and <f> = ir/2. It is seen that the Pareto fronts for <fi = are shifted significantly 
to the left on the time axis (and have correspondingly smaller values of T*) as compared to the fronts for <f> = tt/2. 
This behavior shows that the global phase is an important parameter in TOC of unitary transformations, consistent 
with the findings of Ref. [10]. We also observed for all target gates considered here that the Pareto front for <f> = tt 
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FIG. 6. (Color online) Distance-time Pareto fronts for the QFT and SWAP gates with global phases = and <j> = 7r/2. Each 
front is denoted by a distinct symbol shape and color, as listed in the legend. Fronts for <j> = n and (f> = 37r/2 (not shown here) 
closely match those for (f> = and (f> — n/2, respectively. The significant difference in front locations (and corresponding T* 
values) demonstrates the importance of the gate's global phase in TOC. The black stars denote the average values of T> obtained 
from the searches reported in Table II. 



TABLE I. Critical times T* (4>m) obtained using the PFT method for the selected set of target gates and four global phase values 

m = m7r/2 (m = 0,1,2,3). 



Gate 


m = 


m = 1 


T*{<t> m ) 


m = 2 


m = 3 


CNOT 


4.12 


4.11 




4.15 


4.09 


SWAP 


4.19 


11.80 




4.34 


11.84 


Vswap 


2.26 


9.84 




2.20 


9.82 


QFT 


5.98 


9.86 




5.94 


9.86 


QFT' 


6.19 


9.92 




5.96 


9.90 


CPHASE(tt) 


4.22 


3.98 




4.21 


3.96 


CPHASE(tt/2) 


2.25 


5.98 




2.38 


5.96 



lies very close to the one for = 0, and the Pareto front for (f> = 3ir/2 lies very close to the one for <f> = tt/2. The 
pairwise separation of T*((f> m ) values for global phases with m = {0, 2} and m = {1, 3} is reported in Table I 
for several target transformations. For all gates except CNOT, there is a significant difference in the T* values 
corresponding to the two pairs of <\> m values. 

The dependence of the critical time T* on the gate's global phase <f> has important practical implications for opti- 
mizations in numerical simulations and, potentially, experiments, where implementing the target gate in minimum 
possible time is often desirable. In a common situation where initial control fields are selected randomly, a search 
aimed at minimizing the distance T> for the target gate with a certain value of <f> may, in principle, converge to a con- 
trol solution that enacts the gate with another value of 0. For the target transformation in SU(iV), the discreteness 
of allowed phase values m of Eq. (33) means that the probability of such an event depends on the distance between 
transformations with nearby values of 0. While the fraction of searches that converge to solutions corresponding 
to a "wrong" value of is negligible for two-qubit gates, it will increase for multi-qubit gates as the number N of 
allowed phase values grows (see Sec. IV below for more details). Furthermore, if the phase-independent distance 
Q of Eq. (12) is employed (which is likely in practical situations where the global phase of a gate is undetectable), 
randomly initialized searches will converge to control solutions corresponding to all <p m values (the probability 
for each m should be equal to 1/N, provided that a local search algorithm is employed and initial control fields 
are generated in a truly random manner). Therefore, in order to make all transformations W(<p m ) reachable, the 
control time T must be not smaller than any of the critical times T*(0 r „), i.e., T > max{T*(0 m )}. 

In particular, due to the difference between the T* (0 m ) values, an undesirable situation will arise if T is selected 
such that T*(4> m ) < T < T*(4> m >). Then W(4> m ) will be reachable, while W(4> m ') will not, and therefore local 
searches converging to W(<f> m ) will attain a desired objective value (e.g., Q < 10 -8 ), while local searches moving 
towards W(<f> m >) will be unable to improve the objective value beyond the limit set by the corresponding Pareto 
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front. If such a situation is encountered in a numerical simulation, it may appear as if a fraction of searches are 
"trapped," which seems to contradict the trap-free topology of the quantum control landscape. However, in reality, 
this behavior arises due to the limitation on the control time (i.e., a constraint on a critical control resource) that 
makes one or more of the target transformations unreachable. Obviously, if the control time is increased to satisfy 
the condition T > max{T*(</> m )}, this spurious "trapping" will be completely eliminated. 

We explored this effect by running 500 D-MORPH searches (with random initial fields), aimed at minimizing 
the phase-independent distance Q for QFT and SWAP gates with T = 9 and T = 10, respectively. These searches 
may approach solutions corresponding to any cf> m (0, tt/2, tt, Ztt/2), with optimization trajectories determined by 
initial fields. As expected based on the T* values in Table I, all searches that converged to W{4> = 0) and 
W((f> = tt) attained desired objective values Q < 10~ 8 , while all searches that moved towards W(4> = tt/2) and 
W(4> = 3ir/2) halted well before reaching the target. Table II presents the percentage of searches approaching 
each target transformation W(<f) m ), along with the associated mean value of Q. For the "trapped" cases (occurring 
for <fi = tt/2 and <j> — 3tt/2), the standard deviation around the mean objective value Q is less than 0.5%, which 
indicates that all these runs failed to reach the target due to the same cause. For these searches, we also computed 
and averaged the corresponding T> values, which are shown as black stars on Fig. 6. They lie exactly on the 
respective Pareto fronts, thus demonstrating that the attainable distance values are determined by the limitation on 
the control time. 



TABLE II. Optimization results for 500 D-MORPH runs performed to minimize the phase-independent distance Q of Eq. (12). 
The target gates are Wqft,2 and Inswap with control times T = 9 and T = 10, respectively. The percentage of searches 
that approached each target transformation W((f> m ) = e^ m W is recorded, along with the associated mean value of Q. The 
difference of the reported percentages from 25% (1/JV probability) is mostly due to the fact that initial control fields are not 
completely random, as they are generated using the parameterized form (22). 



Gate 


T 


4>m 


% of searches 


G 


QFT 


9 





28.5 


9.87 x 10" a 






tt/2 


8.6 


4.77 x 10~ :i 






TT 


15.6 


9.91 x 10~ 9 






3tt/2 


47.3 


4.71 x 10 -3 


SWAP 


10 





43.3 


9.95 x 10 _a 






tt/2 


26.4 


2.14 x 10~ 2 






- 


22.9 


9.91 x 10~ 9 






3tt/2 


7.4 


2.13 x 10~ 2 



Further insight into the difference between searches that approach transformations corresponding to <f> m = 
{0, tt} and <p m = {tt/2, 3tt/2} can be gained by examining the path length difference A(s*) — A(s). From 500 D- 
MORPH runs listed in Table II for the gate Wqft,2 and T = 9, we selected 20 typical searches: five for each value 
of 4> m . In Fig. 7, A(s*) — A(s) is plotted as a function of Q for these 20 searches. The optimization trajectories 
corresponding to <f> m = {0, tt} (shown as black dashed lines) successfully converged (Q < 10~ 8 ) to the respective 
transformations and exhibit a behavior similar to that of the trajectories with T > T* in Fig. 2. Specifically, 
A(s*) — A(s) changes with Q according to a power law (which appears as a linear change on the log-log plot) 
over the range 10~ 7 < Q < 10~ 3 . In contrast, the optimization trajectories corresponding to <fi m = {tt/2, 37r/2} 
(shown as red solid lines) halted at a suboptimal objective value Q w 0.0047, near which the "trapping" of these 
searches is manifested by a very large change in A(s) (over ~ 4 orders of magnitude) that produces negligible 
improvement in fidelity. These trajectories also exhibit a significant "flattening" at Q ~ 0.1. The qualitative 
difference between the trajectories corresponding to the converged searches and those destined to be "trapped" due 
to unreachability of their target transformations can be utilized to identify the latter ones at an early stage of the 
D-MORPH optimization. We did not observe any search whose trajectory would indicate that it was attracted to 
an unreachable transformation (corresponding to <fr m = tt/2 or 37r/2) en route to finally converging to a reachable 
one (corresponding to (j> m = or tt). 

We also verified that for sufficiently large T values (i.e., T > 9.86 for the QFT gate and T > 11.84 for the 
SWAP gate), all searches attained Q < 10~ 8 . These results demonstrate that the critical times corresponding 
to all values of the target gate's global phase have to be accounted for in order to ensure that all relevant target 
transformations are reachable at the selected control time. From a practical point of view, due to the extremely 
high search effort in the vicinity of the critical time, it is preferable to use a T value that safely exceeds all T*{<j> m ). 
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FIG. 7. (Color online) Path length difference A(s*) — A(s) along optimization trajectories seeking Wqft,2 with T = 9, as 
a function of the phase-independent distance Q. Black (dashed line) trajectories: searches that attained Q < 10~ 8 as they 
converged to W(<j} m ) with cj> m = {0, 7r}. Red (solid line) trajectories: searches that were "trapped" at Q m 0.0047 as they 
approached W(4> m ) with m = {tt/2, 3n/2}. Trajectories for the converged searches behave similarly to those with T > T* 
in Fig. 2. On the other hand, as the red trajectories approach the "trap" at Q ~ 0.0047, they show a very large change in the 
field structure (A changes by ~ 4 orders of magnitude) with negligible improvement in fidelity. 




C. Dependence of critical time on inter-qubit coupling strength 

For two- and multi-qubit gates, the implementation speed is limited, in most practical situations, by strengths 
of inter-qubit interactions. To explore the relationship between the critical time T* and the coupling strength 
J = jC 1 ' 2 ) for a given target gate W, the PFT procedure was performed with various values of J for an otherwise 
fixed two-qubit Hamiltonian. The resulting distance-time Pareto fronts for systems with 0.2 < J < 3.2 and 
W = Wqft,2 are shown in Fig. 8. As J varies, the front's shape remains approximately the same (on the log-log 
plot), but the entire curve is shifted on the T axis. In Fig. 8, the coupling constant is doubled for each successive 
Pareto front from J = 0.2 (the rightmost curve) to J = 3.2 (the leftmost curve), and the equal spacing between 
the fronts on the logarithmic-scale abscissa suggests that T* has a power law dependence on J with a negative 
exponent. 




T(arb. units) 



FIG. 8. (Color online) Distance-time Pareto fronts for selected values of the inter-qubit coupling strength J. The target gate 
is Wqft,2. Each Pareto front is shown by a different symbol shape and color, with the J value indicated to the right of the 
corresponding curve. The equal spacing between the fronts on the logarithmic-scale abscissa indicates a power law dependence 
of T* on J. 

Figure 9 confirms that T* scales as 1/ J for a significant range of coupling strength values and various target 
gates. Specifically, we considered QFT, QFT', CNOT, SWAP, the identity transformation I4, and a random trans- 
formation W ln d € SU(A^) produced by W Tn d = e%A , where A is a random, traceless Hermitian matrix. For each 
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FIG. 9. (Color online) Critical times T* identified using the PFT method are plotted versus the inter-qubit coupling strength 
J. Points corresponding to different two-qubit target gates are denoted by distinct symbol shapes and colors. Regression lines 
obtained from the least squares fit have slopes in the range from —0.980 to —1.015 and are shown over the ranges of J values 
for which the power-law relation holds. The scaling T* oc 1/J breaks down when J > 0.2 for the identity transformation, 
J > 10 for the CNOT gate, and J > 20 for the other gates. 



of these target gates, the power-law exponent was evaluated as the slope of the corresponding regression line on 
the log-log plot of Fig. 9. These slopes range from —0.980 to —1.015, in excellent agreement with the scaling 
T* oc 1/J. As the regression lines for various choices of W have different intercepts on the ordinate, the prefactor 
in the relationship T = c/J is a function of the gate: c = c(W). The value of the prefactor c for the identity 
transformation is much smaller than for any other gate, while the largest value of c (among the gates considered 
here) is found for the random transformation. 

In previous studies of TOC [2, 3, 9, 10, 33, 34], the scaling T* oc 1/J has been assumed to hold for various 
two- and multi-qubit gates. However, an explicit numerical validation of this property has not been previously 
demonstrated over large ranges of J values. Figure 9 shows that this scaling breaks down when J > O.Olwi 
for the identity transformation, J > 0.5wi for the CNOT gate, and J > uj\ for the other gates (recall that 
u>\ = 20 is the smaller transition frequency of the two qubits). The breakdown of the 1/J scaling for non-identity 
gates happens when inter-qubit and intra-qubit dynamics can occur on roughly the same time scale, i.e., when the 
rate of the quantum information exchange between qubits stops being the "bottleneck" for the speed of the gate 
implementation. 

Distance-time Pareto fronts for strong inter-qubit couplings (J ^> uii) exhibit a distinct structure as compared 
to those for smaller J values. In Fig. 10, we compare the Pareto fronts for (a) weak coupling: J = 0.8 = 0.04a;i, 
and (b) strong coupling: J = 400 = 20wi, obtained for three target gates (CNOT, QFT, and SWAP). To facilitate 
the comparison, the distance T> is plotted versus the reduced time T/T* . When coupling is weak, the fronts for 
different target gates have similar, but distinct shapes. When coupling is strong, the fronts for the QFT and SWAP 
gates exhibit non-monotonic behavior, as the target transformation is unreachable in a range of control times larger 
than T* . Such behavior was also observed in other multiobjective optimization studies [67, 68]. 

According to the data presented in Fig. 9, in the strong coupling limit (J ^> ui), the critical times T* saturate 
at values on the order of ir/uii w 0.157. A possible explanation for this effect is that, in the considered model 
systems [c.f., Eqs. (18) and (20)], one-qubit rotations around the z axis (or, in fact, any arbitrary axis) cannot be 
accomplished in arbitrarily short time because the control fields are polarized only in the x direction. The obtained 
values of T* suggest that a target transformation would become unreachable when the control time is too short to 
perform a phase flip (a ir rotation around the z axis) of the qubit with the smaller transition frequency. A similar 
situation where the critical time is set by the minimum transition frequency was also encountered in state control 
[69]. That work [69] also reported an analytical form of the Pareto front, for a variety of controlled quantum 
systems. The distinct Pareto front shapes shown on Fig. 10 indicate that such a universal analytical form likely 
does not exist for control of unitary transformations. 
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T/T* 

FIG. 10. (Color online) Distance-time Pareto fronts for the target gates Wcnot, VKqft,2, and Wswap (denoted by red circles, 
blue diamonds, and green squares, respectively). The normalized distance T> is plotted versus the reduced time T/T*. (a) The 
fronts for J = 0.8, where V < 10~ 8 is attainable for all values of T > T*. (b) The fronts for J = 400, where the QFT and 
SWAP transformations are unreachable in a range of time values for T > T* . 



IV. CRITICAL TIMES FOR CONTROL OF THREE- AND FOUR-QUBIT GATES 

In this section, we explore TOC of unitary transformations in three- and four-qubit systems, with the primary 
goal of identifying critical time values as the number of qubits increases. In particular, we investigate whether 
the scaling T* oc 1/ J that was demonstrated for a range of J values in the two-qubit case also holds for multi- 
qubit gates. In all simulations reported in this section, we consider the target gate Wqft> „ [Eq. (30)]. In most 
simulations, the inter-qubit coupling strengths in Eq. (18) are selected equal, i.e., j( fc j) = J, Vfc, j. However, 
several cases of unequal couplings are considered as well for three-qubit systems, with the focus on two physically 
realistic scenarios: (i) j( fc ^) values are similar but not equal (e.g., spins on a lattice or in a polyatomic molecule) 
and (ii) J* 1 * 2 ) = j( 2 < 3 ) > J^ 3 ) (e.g., a linear spin chain). The choices of unequal j( fc J) values along with the 
average coupling strength J = [J^ 1,2 ' + j' 1 ' 3 ' + j' 2 ' 3 )]/3 are displayed in Table III. 

Critical times T* obtained from PFT trajectories are plotted in Fig. 1 1 versus the coupling strength J (for equal 
couplings) or versus the average coupling strength J (for unequal couplings). For two-, three-, and four-qubit 
systems, the dependence of logT* on log J, in the range of coupling strengths J < uj\, is evaluated using the 
least squares fit. The slopes of the resulting regression lines on the log-log plot of Fig. 1 1 are —1.015 ± 0.002 for 
n = 2, —1.00 ± 0.02 for n — 3, and —0.99 ± 0.07 for n = 4. Thus, in the range of validity, these numerical data 
agree with the scaling T* oc 1/ J very well for each n. Also, the T* values for three-qubit systems with unequal 
couplings follow (for J < u)\) the same scaling with respect to J (i.e., T* oc 1/J), and, overall, these points fit 
very well with those obtained for equal couplings. Thus, critical times for systems with arbitrary j( k >j> values may 
be reasonably well estimated from equal-coupling results with J = J. 

The data obtained for both two- and three-qubit systems demonstrate that the scaling T* oc 1/ J breaks down for 
J > = 20. As J increases beyond u>\, the T* values for two- and three-qubit systems gradually converge and 
trend towards tt/lji 0.157. As noted in Sec. IIIC above, this value coincides with the minimum time required 
to perform a phase flip of the qubit with transition frequency u>\. At J = 400 (the strongest coupling that we 
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TABLE III. Values of inter-qubit coupling strengths j( k 'i> employed in simulations with three-qubit systems. The average 
coupling strength J is also shown. The corresponding T* values are plotted versus J in Fig. 11. 
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FIG. 11. Critical times T* plotted versus the inter-qubit coupling strength J, for n-qubits systems (?i = 2,3,4). Different 
system sizes are denoted by distinct symbol shapes and colors. The target gate is Wqft'.ji- Regression lines obtained from 
the least squares fit are shown for systems with equal couplings, and their slopes range from —0.99 to —1.015. The scaling 
T* oc 1/J breaks down when J > 20. For three-qubit systems with unequal couplings, T* values are plotted versus the 
average coupling strength J (magenta diamonds), and these points fit very well with those obtained for equal couplings (blue 
circles). 



considered), this critical time value is already reached for n — 2 (T* « 0.157) and is very closely approached 
for n = 3 (T* « 0.162). The saturation of T* at nearly the same value for both two- and three-qubit systems 
further suggests that, in the strong coupling limit, the critical time is determined by the minimum time required to 
accomplish one-qubit rotations. 

When the system size increases, the difference between nearby allowed values of the global phase decreases 
as \4> m — 4> m ±i\ = 2-k/N, and the normalized distance between the corresponding transformations decreases 
as T>[W(4> m ), W(4> m ±i)] = sin 2 (7r/7V). Therefore, when the phase-dependent distance T> is employed as the 
control objective for three- and four-qubit systems, a small fraction of randomly initialized searches may converge 
to transformations that differ from the target one by a global phase. Indeed, we encountered such a situation in our 
simulations, when a small number of D-MORPH searches aimed at reaching the target gate Wqft' 3 or Wqft',4 
converged instead to gates e ±47r / 4 l4 y QFT'.3 or e ±47r / 8 M / QFT',4> respectively. As discussed in Sec. IIIB, when such 
convergences to gates with nearby values of <j> happen, a sufficiently large control time T > max{T*(0 m )} is 
required to ensure that all transformations W((f> m ) are reachable. If this condition on T is not satisfied, at least 
some local searches will be unable to attain a desired objective value due to the character of the distance-time 
Pareto front below the critical time. Such a situation where a small fraction of optimization runs employing the 
BFGS method were seemingly "trapped" was reported in Ref. [58] for a three-qubit spin-chain system and the 
QFT target gate. This effect is fully explained by the fact that the optimizations in [58] used the T value such that 
T*(4> m ) < T < T*(cf> m i) [70], thus making transformations corresponding to some phase values unreachable. 
We reproduced the results of Ref. [58] and verified that this spurious "trapping" is completely eliminated when the 
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control time T is made sufficiently large. 



V. CONCLUSIONS 



This work examined TOC of quantum unitary transformations through the exploration of the Pareto front that 
quantifies the trade-off between the goals of minimizing the distance to the target gate and the control time. The 
PFT algorithm was introduced to (1) identify the critical time T* below which the target transformation is not 
reachable and (2) move along the Pareto front to find families of optimal control fields that minimize the distance 
(or, equivalently, maximize the fidelity) at various values of T. Our results suggest that a distinct Pareto front exists 
for each selection of the control system and the target gate. A practically relevant feature observed for many gates 
is the strong dependence of the critical time on the global phase of the target transformation. We also examined 
the dependence of T* on the inter-qubit coupling strength J and confirmed the universal scaling T* oc 1/J (or, 
T* oc 1/ J for unequal couplings), consistent with expectations in the literature [2, 3, 9, 10, 33, 34]. However, we 
found that, while this scaling holds for a wide range of J values, it breaks down when J > wi, where cji is the 
smallest of the qubit transition frequencies in the system. Thus the critical time cannot be made arbitrarily small by 
increasing the coupling strength, and the ultimate limit on the T* value is set by the smallest transition frequency 
oji which, for the considered model systems, is the minimum speed of one-qubit rotations. 

One of the goals of this work was to relate the obtained TOC results to properties of the optimal quantum control 
landscape. In particular, we observed that, for a given quantum system and a target transformation, different 
randomly initialized PFT trajectories produce essentially the same Pareto front. This result suggests that (1) as 
long as T > T*, the favorable landscape topology characteristic for unconstrained control fields remains intact, 
and (2) for T < T*, the limitation on the control time significantly affects the value of the landscape's global 
optimum (so that unit fidelity/zero distance is no longer attainable), but does not lead to landscape fragmentation 
and formation of local traps (provided that we distinguish between landscapes for target transformations that differ 
by the global phase). We also found that the optimization search effort rises superexponentially as T decreases 
and approaches T*, with corresponding changes observed in metrics quantifying the local structure of the control 
landscape. In particular, "flattening" of the landscape near the optimum correlates remarkably well with the search 
effort growth as the control time decreases. 

The structure of the distance-time Pareto fronts identified in this work has important implications for experimen- 
tal implementation of high-fidelity quantum gates in short times. Taking into account the very steep slope of the 
Pareto front at T < T* and realistic experimental uncertainties in the gate control time, operating at a T value that 
safely exceeds the critical times corresponding to all allowed values of the target gate's global phase is desirable 
to avoid a significant loss of fidelity. Furthermore, values of T very close to T* should be avoided in order to keep 
the optimization search effort at a reasonable level. 

In summary, this work extensively investigated the effects of limiting the time as a resource for optimal control 
of quantum unitary transformations. While our study employed a particular physical model, the above results 
are expected to be qualitatively applicable to other coupled-spin systems. In addition to the goal of generating 
target unitary transformations, the methods presented here can be applied to other objectives in quantum control, 
including state preparation and optimization of observable expectation values. 
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